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abstract 

The Galerkin weighted residual technique using lin- 
, r i angular weight functions is employed to develop 
finite difference formulae in Cartesian coordinates for 
hP LaolMiM operator on isolated unstructured trian- 
gular grids The weighted residual coefficients asso- 
rted with the weak formulation of the Laplacian 
operator along with linear combinations of the 5 esi *|® 
Si reused to develop the algorithm. The algo- 
was tested for a wide variety of uns ructured 
meshes and found to give satisfactory results. 


Subscripts 

e+ index, Eq. (8) 

fd finite difference or Taylor series 

j element index, Eq. (10) and Mg. 1 

k element index, Eq. (10) and Mg. 1 

M total number of elements 

a global node index, Eqs . (6) and (7), Fig. I 


nomenclature 


nm index, Eq. (9) 


A area 

A e element area, see Fig. 1 

Aj total cell area, see Fig. l 

A V Voronoi (finite difference) area 

B bias constant, Eq. (13) 

f linear combination factor, Eq . (22) 

h coordinate spacing distance 

M number of elements in cell 

N linear interpolation function 

W weight function, see Eq. (2) 

x axial distance coordinate 

y transverse distance coordinate 

a weighted residual coefficient, Eq. (2) 

[3 property constant, Eq. (I) 

4 ) dependent variable, Eq. (1) 


□ central node index 

p number of outer nodal points, see Fig. 1 
a index, Eq. (11) 

index, Eq. (5) 

Supersc r ipts 

- numerical estimate 

average value , Eq . (. B3 ) 

(e) element value 

n global element index, Eq. (3) 

run global element index, Eqs. (3) and (8) 


INTRODUCTION 

Linear triangular finite element theories 

(Baumeister and Horowitz, 1984) ^MpoMed^as U Y 
applied to the Fourier transformed linearized ga 
dynamics equations (-steady" state equat tons) in t the 
study of noise propagation in variable area ducts 



the acoustic design of turbofan engine nacelles c im ; 

Hor/h transi ® nt S0 ^ ut i°ns to the acoustic wave equa- 
tion have used conventional finite differon^ ^ 

and elliptic mapping methods (Raad and White 6 iq86Mc> S 

' f i 1 " very u '*' «■<*•! «;r 

£«;» ;; the 
rans tent technique employing elliptic mapping has been 

unsuccessful in pred.cting the acoustic far ffeld radi 
anon pattern from aircraft inlets. To overcome this 
problem the present study is concerned with adap ing 
n sly developed 'steady" state triangular fin, 
ement mesh generators (Baumeister and Horowitz 1984) 
for use transient variable geometry acousUcs o 
SlEE; the f0rm unstructured ffiiS 

In developing unstructured finite difference eoua 
ons, Jameson ( 1987 ) ( Eq . (4. 10)) has applied the Q 

time d Hpn T 8 ^ residual Calerkin method to obtain a 
time dependent discretization of the Euler equations 
for an unstructured triangular mesh. Erlebacher (1985) 

to establish a difference algorithm for the LaDlarian 

spec? i vp F r!.' I 1 Unf ° r ^na tely , from a pointwise per- 
specfive, Erlebacher has shown that the difference 

on?v r ?prr> S ’a 3 ' result from a y ariat ional approach are 
, l ^ e [° order accurate (Erlebacher, 1985 . p 47) f or 
the Laplacian operator on nonuniform unstructured grids 
as well as many uniform grids. Baumeister (1988) 8 
° b * a ‘ ned f lm ^ lar results using the Calerkin weighted 
residual techniques. The Galerkin algorithms n 

= f ° r the bap,ac ' a " operator w«e con- 

sidered valid for a general unstructured grid when used 
only in con)unct ion with a complete global grid system 
as commonly employed in conventional f ini tenement 
HempnJ 5 ' , Baumeistar U 988 ) also showed that finite 
are related Jv"’* 6 differ f ce ( Ta Vlor senes) algorithms 

• Igor? tali ft ??r«“,;E' a e £,£*! ,he n " ,te 
global average. V senes ap P roach on a 

At present, the difference expression for the 

in i local°vr r H t0r | baS ^ ° n variat,onai theory is valid 
when the irL ?,,? 7 ^ J es ' ricted conditions, namely 

generate finitP H ff S satlsfled (Baumeister, 1988 ). To 

res,d^i n equatiois h Unear COmblnations of ‘he weighted 


GALERKIN WEICHTED RESIDUAL FORMULATION 

sti.dv T !!f bapla ^ ian operator plays a major role in the 
Consider the^oUn 3 ^ e 1 eG 1 romagne t i c wave propagation, 
operator followin S expression for the Laplacian 




( 1 ) 


is 


where <fr is a scalar or potential quantity and 8 

a variable property coefficient that depends on its 

spat. al orientation. To obtain the finite difference 

la X be r L S d Sl T a °r ^ L / P A aci ™ operator at some pom' 

nec'ing the P ^ !i I m* 6 " 11 15 defined b y con 
a , no< ^ a (grid) points surrounding the 

«!? TlLl * s sh r in Fis - >• tAS?,. 

rerence cell of total area A T is then divided into 

?irid) do mis" 8 ? ' ar h 3S Ae staJted out ^ nodal 
(grid) points P as shown in Fig. 1. The number of 


< X 1'V1> 



FIGURE 1. - FINITE DIFFERENCE CELL 
CONFIGURATION. 

connected areas (called an element of the cell) needed 

labeTed n % * d ‘ fferenCe ^ uatlon node o is 

Using the weak formulation of the method of 
weighted residuals (Lapidus and Pinder, 1982 p 443) 

ihLin P ! aCla g ope ^ ator ln Et l- (1) can be integrated to 
obtain (see Appendix A for brief derivation) 


at- 




3V<t>j dA = a Q^g 


¥1 


a 4> 
n y n 


S*-) 


a 6 
P^P 


( 2 ) 


Using the Galerkin approach, the spatial weisht W i « 
JfST'TS by ,h, standard line,? I„?„^!fv[on tap- 

?s the e hat T r the P° tentlal * indicates that 
t , 1 ! * approximate numerical solution for the poten- 

- 1 : *n represents the potential at any general 

noda point n. The value of the general L coeff,- 

U en iVanH n ?h ° nly 00 the loCation of the 8 r ‘d Points 
fo?m >n ° he pro P erty coefficient 0 and are of the 

* (x n - V»V - »„>] - [6 ( " l /aj'")][(y o _ y nJ 
X| V - y„J * (V - *„)(«„ - x„i] m 

“° ‘ ' a?i [ Bte>/ - ,a<e) ] [(y e - y e ,) 2 » (x et - x e ) 2 ] (4) 

The indices are defined as 


, p ~ 9 


* p or 


, p “ 1 

a = p 

(5) 

n- = n- 

1 

+ Ps n,l 


(6) 

n+ * n + 

1 

+ P5 n,P 


( 7 ) 

e+ = e + 

1 

+ Ps e,P 


(8) 

nm = n - 

1 

+ Ms n. 1 


(9) 


p lt 110 1 (Eq nt (3 r i a i)) iVen by Lap ‘ dus ’ and Pinder - 1982, 
i(e) 


bJt x o (y j - yk> + x j(yk - y 0 ) + x k ( y 0 


•] 


( 10 ) 


LAPLAC I AN OPERATOR 

The present paper will establish a relationship 
between the Laplacian operator and the * coeffici- 
ents from the method of weighted residuals as given by 
Eqs. (3) and (4). An expression for the Laplacian ot 

the form 

- f ( a t . ' • • • ) J 

is desired. Using the result from Appendix B for a 
linear shape function, the Laplacian can be expressed 

as 


where At is the total area of the cell. 

For the six point hexagon difference grid systems 
as shown in Figs. 2 and 3. the difference equations tor 
the Laoiacian as calculated from Eq. U2) are in agree- 
ment with the standard Taylor senes finite difference 


operators that appear in the literature iRektorys. 

1969. p. 1114). However, this will not generally oe 
the case, as discussed earlier m relation to Eribacher 
and Baumeister s works and also discussed in Lapidus 
and Pinder 11962). p. 105. For uni form gr ids . of the 
type displayed in Figs. 4 to 5. the finite difference 
and finite element expressions for the Laplacian opera- 
tor will differ with the conventional Taylor senes 
difference equation. Only when the bias constant B 
defined by the following area rule is identical to zero 
will the finite element algorithm and Taylor series 
approximations be equal. 



Ay is the Voronoi neighborhood {Eribacher, 1985) or 
the area normally employed in the finite difference 
analysis and is the area enclosed by the dashed line 
in Figs. 2 to 5. Baumeister (1988) presents a simple 
algorithm to determine the Voronoi {finite difference) 
area for an unstructured grid based on the a coeffi- 
cients in Eqs. (3) and (4). „ 

[f B is not zero in Eq , (131. the coefficients 
in front of the weighted residual expression for the 



FIGURE 2. - UNIFORM 6 NODE HEXAGON GRID. G R [D . 



FIGURE 4. - GLOBAL CONVERGENCE EXAMPLE WITH RECTANGULAR ELEMENTS. 
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FIGURE 5. 
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2 x ] ♦ { -0. •'!) x B - 0 


B ; O.A A r - 5 A y 1 
GLOBAL CONVERGENCE WITH TRIANGULAR GRIDS. 


Lajjlacian operator will differ from the expression for 
the Taylor senes expression by a factor of 1 plus fl. 
it fl is zero, a single pattern of uniform elements 
can normally be placed in a regular domain. For the 
uniform grid systems in Figs. 2 to 5, the area rule 
leads to the fol owing simple relationship between the 

?Fn gh ion eS ‘a Ua u expression for the Laplacian operator 
(tq. (U)) and the standard Taylor series (finite dif- 
ference) representation 


7-$ 


3A„ 

v 


| 


fd 


(14) 


Figure 6 displays the difference equations devel- 
oped from Eq . (12) for a nonurn form six-node cell. 4s 
discussed previously, the Laplacian operator shown in 
rig. 6 is valid when taken in conjunction with the 
other grid systems that surround it in a given finite 
element type domain. However, since the bias of 
-O.U 6 ? is different from zero, the difference equation 
displayed in Fig. 6 will not be valid for a local iso- 
lated cell. For nonuniform grids, the fol lowing test 
functions will be used to check the validity of the 

difference formula on a local unstructured nonuniform 
grid. 


For the conventional four point square difference arid 
as shown in the upper right of Fig. 4, the difference 
equation for the Laplacian is different by a factor of 
3/. from the conventional Taylor series representation 
ot the Laplacian operator (Rektorys, 1969, p. 1114) as 
predicted by Eq , ( 14 ). Similarly, the eight-node sys- 
tem shown in the upper right in Fig. 4 differs by a 
factor ot 3/4. However, since the method of weighted 
residual (as applied in finite element theory) is guar- 

an a e a d ,i!° conver 8 e ; >hese are acceptable formulae pro- 
vided they are applied to a consistent global mesh and 
t a single arbitrary cell. The purpose of the 
present paper is to adapt Eq. ( 12 ) so that it can be 
apphed on a local grid independent of the surrounding 

Although the weighted residual difference equation 
tor the Laplacian operator in Fig. 4 differ in each 
cell from the standard Taylor senes difference equa- 

wp?oh t pH° m a . global point of view the Taylor senes and 
weighted residual equations will average out to the 

wfl? k! nn a'"", T ' s not 2ero ' the global domain 
h p fl ' ledwtth different element types, as shown 
in the combined mesh of Fig, 4 . In the case shown in 
t* l g. 4, the area average of the eight-node element com- 
bined with the four-node element will be such that the 
Taylor series will be valid. As shown in Fig. 4. the 
average global value of B is equal to zero. IB is a 
measure of the difference between the finite element 
and Taylor Senes as just defined.) A similar situa- 
tion occurs for triangular grids as shown in Fig. 5 


xy 


y“ 


( 15 ) 


The difference equation in Fig. 6 satisfies the 1 x 
and y test functions; however, it does not satisfy 
the xy x-, and the y- test functions. It does how- 
ever satisfy an x- * y- test function which leads to 
a simple formula for the Voronoi area in terms of the 
<x coefficients (Baumeister, 1988, Eq, (10)). 

Equation (12) will now be adapted to produce dif- 
terence equations valid on a local grid. That is. a 
difference equation which satisfies the six constraints 
given in Eq. (15) . 

MODIFIED LAPLACIAN OPERATORS 

The method of linear combinations will be used in 
the next section of this paper to adapt the method of 
residuals to obtain local convergence for uniform grids 
and unstructured grids of general shape. For conven- ‘ 
lence. a modified version of Eq. (12) will be used in 
the combination process. 

As shown by the area rule, if the bias constant is 
equal to zero (B = 0) , then Eq . (13) becomes 


A T 


r = A v 


( 16 ) 


Equation (16) combined with Eq . (12) yields 


„2 Vo * “1*1 * V? * • • • 

♦ — (B - 0) (17) 


4 


9 


7^ 

v d) 



FIGURE 


6 . - NONUNIFORM 6 NODE GRID SYSTEM. 


where A u is the Voronoi area of the cell and the a 
values are given by Eqs . (3) and (4). Surprising y, 

Eq (l?) was checked against a number of uniform grids 
such as Shown in Figs. 4 and 5 where B was not equal 
to zero and was found to be identical to the Taylor 
series representation which satisfies the conditions 

given by Eq. (15), *. 

In the linear combination process to follow, aroi 
trary constants will be assigned to the Laplacian oper- 
ators. Therefore, either Eq. (12) or (17) could be 
used in the weighting process for unstructured grids 
since they only differ by a constant for a given dif- 
ference cell. However, since Eq. (17) fortuitous ly 
seems to work for all uniform grids system, Eq. (17) 
will be employed in the superposition process. 

LINEAR COMBINATION FOR UNSTRUCTURED GRIDS 

The method of linear combinations has been 
employed in the past to obtain higher order difference 
equation for finite difference grid systems. For 
example, Salvadori and Baron (1962), p. 235, used lin- 
ear combinations of a four-node system to obtain the 
difference equation for the eight-node system shown in 
Pig. 7. In this case the difference equation for the 
eight-node system can be written as 


7“4> = - f : 


( 18 ) 


Equation (18) automatically satisfies the first tour 
conditions of Eq. (15) since operators 1 and -individ- 
ual I v -neet these conditions ot a zero valued iap acian. 
Substituting either the x- or the y- term into 
Eq. (18) requires a Laplacian value ot . or .nat 

1 = fi -f 2 U3) 

Thus, there are an infinite number of combinations of 
difference equations that can represent the eight -node 
grid system in Fig. 7. For example if U ■ J“W» 

= 1/3 the difference equation for the Laplacian 
becomes (Rektorys. 1969. p. nil) 


7“(t> 


4($, + 4>, i- 41, + 4> 7 ) * (4> 2 ► 4»^ + + V ~ 30<t 0 

• ~ “ T 

6h“ 


or if fi = 1/3 and f 2 - 2/3 the difference equati 
for the Laplacian takes on the form 


( 20 ) 

on 


(3) 


CM) 


( 2 ) 


(M) 


(31 


( 2 ) 



2 d) 1 + <t> 3 + <d 5 f <t> 7 - l % 

# 



(5> 


► < 

.. c 

►——41 

■v(O) 

o 

( i < 

V “ " 

1 n 


( 1 ) 


( 6 ) 


a) 

v 2 o 


( 8 ) 


FIGURE 7. - LINEAR COMBINATION OF OPERATORS USING FINITE DIFFERENCE THEORY. 
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*6 f *7 


^8 - 8 * 0 


This is the ident ical expression obtained from finite 

hTV^Lo US l " 8 blllnear we ‘§ht functions (Lapidus 
and Pinder 1982 p. 104). Although both equations are 
equally valid, the round-off errors associated with 
each equation are different. An evaluation of the 
round otf error can be found in reference (Salvadori 
and Baron, 1962, p. 236). 

will exajI, P 1 ® 1 th e linear combination method 

If! u 6 !° d ? er " nne the Laplacian operator for the 

i h 4 nd r Ce i.M° Wn u ln Fig ' 8 ' Thls raain ce| l is 
6d 1 w hile the f ? and f 3 cel Is are subd.vi- 
sions of the main cell using only a portion of the grid 
points to define the cell. The Laplacian operator for 
cell l can be expressed as the linear combination of 
hese cells in the following form with f's representing 
the unknown linear combination coefficients: 

7 * ’ f l 7 ? + f 2 7 2* + f 3 7 3* (22 ) 

Equation (22) satisfies identically the first three 
conditions in Eq. (15) for any values of the f con- 

Iatp! S th‘ nCe Eq n- (17) satisfies 'hese conditions. The 
later three conditions in Eq. (15) will be satisfied 
by a simultaneous solution of the following three con- 
straint equations with unknown f, , f-, and f? 
coefficients: J 

(<f> = xy constraint) 

■ ° < 23 > 

(<J> = x- constraint) 


W| , + wi * f 3 7^i 
W - - W 3 3 I 


( 4> = y 2 constraint ) 

, + f 2 7 2»| , + f 3 7 ?| 


The Miown values of the Laplacian terms n Eqs r’M 
to 1251 1 are easily determined from Eq . . in with 7 he 
i coerticients evaluated from the known position of 
the grid points is required m Eos. * ,| ) and | -i | "he 
values used to evaluate the Laplacian terms' V ' 

tac,!r,« •>« ai 1.1.1 

in JiTfS/r Eqs: 1231 “ 1251 =“ » »..«•« 


7 ?j. 

7 ;o 


7 > 


f=xy 


$=xy 

0 

i-xy 

7 ! 9 ! , 
f>=x- 

7 ? 


0 

7 ? 

. . 


7 2*f 


7 ?| 



S ,0n , (36 ! is a ver V simple system easily solved by 
third-order determinants for the values of fj. f, , a „d 

The Laplacian difference equation for the six-node 
grid in Fig. 8 is shown by the formula in the lower 
portion or the figure. This equation meets all the 

f p n" 1 ; 1 ^ g>ven by Eq . (15). Other subdivisions 
of cell 1 could have been performed to obtain a dif- 

ferenorm of the difference equation. Similar to 

Eqs. („o) and (21) it is expected that these equations 
would also have different values of the round-off 
error . 

Similar results for the five-node cell system are 
shown in Fig. 9 . An attempt to adapt the approach to 
he four-node system in Fig. io proved unsuccessful. 
Apparently, the matrix equations are linearly dependent 
in this case and a solution is not possible. 

HIGHER ACCURACY 

Increasing the number of grid points in the dif- 
ference operator generally improves the accuracy of 
finite difference algorithms. A similar improvement in 

the re^ U . aCy I! 0U d u® obtain f rom E Q • (26) by extending 
the relationships shown in Eq. ( 15 ) to include higher 

order terms in the Pascal polynomial triangle. The 
Pascal s triangle for complete polynomials in terms of 
two variables x and y can be written as: 



i. C 

V 2 <t> = 0.264000, + 0.0349902 * OJWS^ ♦ 1.39700,, + O.SSOO^ + 0.805390 6 - 0.160SO 0 
FIGURE 8. - SIX NODE LOCAL LAPLACIAN OPERATOR. 
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( 3 ) 



( 1 ) 



( 2 ) 


f v ^0 


W* 

= -0.053815®, * 0.8656® 2 ♦ 0.39222® 3 - 1.79210,, - 1-29840^ 
FIGURE 9. - FIVE NODE LOCAL L API. AC IAN OPERATOR. 


f 3 V 3 2 ® 


4 . 3 145011 



xy 


x 3 y x 2 y 2 X y 3 y 4 


(27) 


If a sufficient number of external nodes exist, 

Eq (26) could be extended to include the cubic or ev 
the quadratic terms in Eq. (27). This should improve 
the accuracy to the difference algorithm by decreasing 
the truncation error . 

CONCLUSIONS 

The Galerkin weighted residual technique using 
linear triangular weight functions and the method of 
linear combinations are employed to develop finite di 
ference formulae in Cartesian coordinates for the 
Laplacian operator on local unstructured triangular 
grids The algorithm was tested for a wide variety of 
unstructured meshes and found to give satisfactory 
results . 


APPENDIX A - DERIVATION GALERKIN WEIGHTED RESIDUAL 
EQUATIONS 

Let 5 be an approximate weighted residual solu- 
tion for the potential * defined by Laplace s equa- 
tion Next, the total cell area A-p is divided into 
non. e.xi snatial dependence 

elements as shown in Fig. l. ne spaii a i y 

» in each element can be approximated by 


M 

of 


i (e) (x.y) - £ N| e) 
i = l 


(x.y 


,^(e) , | N ‘ e 'j{* (e) } (AD 


3) 


(in Fig. 1:0-1. j - 2. 

The standard linear local element shape function^) 
is employed here. These common linear pyramid like 
interpolation functions S are presented n most 
c : n i ip p I pnipn t texts ( Lap i dus and Pinoer, 198—. p. 

Eq (3 1.4)). The superscript (e) is used to designate 
the element value. Likewise the superscript e) on 
the it designates the i nodes associated with a 
particular element. 


7 


The '.veak formulation of the weighted residual 
approach is used (Laptdus and Pinder. 1982 p 443 
integration by parts plus the divergence theorem) along 
w'th ,ne Gaierkm approximation to obtain the standard 
te.x book weignted residual approximation for the 
Lap I ac 1 an : 


M 

£ 

e - 1 


le ) 


1 

■j j 


i -7N 


(e) 


A 


7 N 




dA = 0 


(A2) 


Equation) A 2 ) can now be easily evaluated using 
conventional finite element theory to obtain the dif- 
ference coefficients listed in Eq . (2) in the body of 
this report for the central cell node labeled o. 

APPENDIX B - AVERAGE OPERATORS 


form 


Consider the partial differential equation of the 


7> 


S<t> 
Sx 1 


(Bl) 


Employing the Galerkin weight residual solution of 
<4- ( 1) tor the four-node cell in the upper right of 
rig. 4 yields 


[atjiflj +■ ayjc, e 


*4*4 


Vo J 


= 0 


The coefficient; associated for the Laplacian have been 
lett in general form while the coefficient of the first 
derivative have been written for the specific grid 
associated with the four-node grid of Fig 4 

it is desired to put Eq . (B2) in a more ' fami 1 iar 
form which coincides with the standard Taylor Senes 

;^ , .! e J 1 L fe T Ce a PP roach - Th e average difference 
alue of the first derivative using the method of 
residual can be written as (Baumeister , 1988 
Appendix B) 


frJ J, N o dA 

T 


i Utj 


< v /3 . 


(B3) 


!“ iVrfn 56 ' the lnteg f al of the wei « h ‘ \ over the 
"f? I? equal to A t /3 , which was determined 

nin-a h a n f ! he standard af ea integration formula 
(Lapidus and Pinder, 1982, p. 112 ) or as shown expli- 

A F Iaire 0985), p. 45). Therefore, in 
Eq. division by the area factor At/3 will con- 

ver the weighted residual coefficients into the 
Taylor senes difference approximations for the first 
derivative. Furthermore, for the specific four-node 
grid shown in Fig. 4. A T corresponds to 2h~, thus 
Eq. ( B 2 ) becomes 


1*1®! 


* 


a 3*3 


V 4 




fp, 


2 h 


^ 0 
J 


I B4 | 


The right -hand term is immediately recognized as 'he 
Taylor series difference approximation for the ff-st 
derivative term in Eq. (Bl) for the four-node grid svs- 
tem shown ;n F, g^ 4 Likewise, the left-hand term can 
cia ed e w, w ih 'n Lap,acian dlff «rence equation asso- 

for Eo -V F 5 ,' l erm f in Ed - |Bl > » ld is the basis 
tor Eq . ( i_ ) m the body ot this report. \ similar 

denvat ion would apply to more general linear second- 
order partial differential such as the Helmholtz 
equation. 
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